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1 .   INTRODUCTION 

Pressure  vessels  and  naval  structures  are  obtained  by  specific  assemblage  of 
plate  and  shell  panels.   The  knowledge  of  the  behavior  of  these  structural 
components  is  essential  if  optimum  design  and  integrity  of  the  overall 
structure  for  a  given  set  of  parameters  are  sought.  [1],  [2] 

This  report  deals  with  the  nonlinear  analysis  of  arbitrary  thin  shell 
structures  subjected  to  static  loads.   The  nonlinear  analysis  includes  pre  and 
post-buckling  behavior  for  any  degree  of  nonlinearity  due  to  large 
displacements  and  large  rotations,  but  small  strains. 

Despite  the  important  research  and  development  efforts  made  since  the 
beginning  of  the  finite  element  method  era,  the  analysis  of  shell  structures 
is  still  an  open  active  research  subject.   The  following  questions  are  still 
actively  investigated: 

-  how  to  approximate  the  real  three  dimensional  problem? 

-  what  type  of  finite  element  discretization  is  most  appropriate? 

-  how  to  solve  accurately  and  efficiently  the  nonlinear  equations  in 
various  situations  of  pre  and  post-buckling? 

-  how  practical  and  general  is  the  computer  code  that  aimed  to  solve  the 
problem  and  what  are  its  computer  resource  needs? 

The  object  of  our  report  is  to  present  a  formulation  which  includes  some 
recent  developments  on  nonlinear  continuum  mechanics,  plate  and  shell  finite 
elements,  automatic  solution  and  strategies  for  nonlinear  equations  and  to 
present  the  possibilities  of  a  computer  code  that  is  adapted  for  mini  and 
micro-computers  to  solve  moderately  small  size  shell  problems. 


The  nonlinear  formulations  considered  are  a  Total  and  an  Updated  Lagrangian 
formulations  [3],  [4],  [5]  combined  with  flat  simple  triangular  elements 
having  only  3  nodes  and  the  6  engineering  degrees  of  freedom  at  the  nodes. 
The  shell  finite  element  is  obtained  from  the  superposition  of  the  CST  and  the 
DKT  plate  bending  element  known  to  be  very  efficient,  reliable  and  effective 
for  all  thin  plate  bending  analysis.  [6],  [7],  [8],  [9],  [10],  [11],   The 
nonlinear  equations  are  solved  using  various  methods  and  strategies  based  on 
the  full  or  modified  Newton-Raphson  method  to  deal  with  the  automatic 
determination  of  the  pre  and  post- buckling  load  displacement  curves.   Three 
basic  strategies  are  considered:   the  load  incrementation,  the  displacement 
control  method  [12],  [13],  and  the  arc-length  control  method  [14],  [15],  [16], 
[17],  [18].   The  FORTRAN  77  routines  dealing  with  the  shell  element  and  the 
nonlinear  solution  procedures  are  compatible  with  the  documented  computer  code 
MEF  presented  in  detail   in  [19].   The  numerical  examples  presented  in  this 
report  have  been  obtained  using  a  VAX  11/750,  a  VAX  11/780  and  an 
APOLLO/DN300. 

The  Lagrangian  Formulations  (TLF  and  ULF)  considered  in  this  report  are 
discussed  in  chapter  1.   The  DKT18  triangular  shell  element  is  described  in 
chapter  2.   The  solution  strategy  to  deal  with  the  automatic  determination  of 
the  load  deflection  curves  is  presented  in  chapter  3.   The  numerical  results 
are  discussed  in  chapter  4.   They  deal  with  nonlinear  behavior,  buckling, 
post-buckling  and  large  rotations  of  elastic  shells  subjected  to  one  variable 
load  parameter. 


2.    THE  LAGRANGIAN  FORMULATIONS  FOR  NONLINEAR  SHELL  ANALYSIS 

2.1   Different  configurations  of  a  shell  in  space. 

We  consider  a  shell  structure  with  a  fixed  orthogonal  coordinate  system 
XYZ  (Fig.  1): 

-  T  refers  to  the  undeformed  (initial)  configuration 

-  ^T  refers  to  a  deformed  (intermediate)  configuration  in  equilibrium 
under  a  given  set  of  loads  *f 

-  T  refers  to  a  deformed  (final)  configuration  in  equilibrium  under  a 
given  set  of  loads  f 

The  purpose  of  the  study  is  to  describe  as  precisely  as  possible  all 
deformed  shell  configurations  like  *T  and  V   for  given  sets  of  loading, 
prescribed  displacements,  boundary  conditions,....   The  description 
includes  def ormational  aspects  (displacements  for  °  T,    rotations,  strains) 
and  mechanical  aspects  (true  stresses  at  the  material  points). 

Two  Lagrangian  formulations  are  considered  in  this  report: 

The  Total  Lagrangian  Formulation  (TLF) .   In  this  case  all  quantities 
(displacements,  strains  and  stresses)  in  the  computational  process 
are  related  to  the  undeformed  initial  configuration.   The 
intermediate  configurations  ^T  can  then  be  interpreted  as  final 
configurations  for  different  sets  of  loads  or  prescribed 
displacements.   To  obtain  the  exact  T   solution  for  thin  shells  under 
Kirchhoff's  assumption,  the  nonlinear  Green-Lagrange  strain  - 
displacement  expressions  must  be  complete  (no  terms  neglected). 
These  complete  expressions  are  very  complicated  and  imply  second 
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Figure  2  Coordinate  Directions  x,y,z 


derivatives  of  all  the  components  of  the  displacement  vector  u. 
Therefore,  this  approach  has  not  been  retained  for  practical 
nonlinear  analysis  unless  further  assumptions  (like  the  moderate 
rotation  hypothesis)  are  made. 

The  Updated  Lagrangian  Formulation  (ULF) .   In  this  case  an 
intermediate  configuration  ^T  is  used  as  a  reference  configuration 
to  obtain  the  (final)  configuration  r  for  a  given  set  of  parameters 
(loads,  ...).   The  *T  configuration  is  supposed  to  be  known  (  *r  is 
in  fact  a  previous  T   configuration).   That  is,  *F  can  be  fully 
described  both  from  the  point  of  view  of  geometry  and  of  mechanics 
(internal  stress  field).   There  is  theoretically  no  difference 
between  TLF  and  ULF  (they  both  want  to  solve  the  same  equilibrium 
problem)  that  is  to  find  T.   But  practically  we  can  take  advantage 
of  the  fact  that  ^T  is  known  and  that  we  want  to  obtain  a 
configuration  T   "not  too  far"  from  T.   In  ULF,  T  is  a  neighboring 
configuration  of  ^T.   Hence,  we  can  consider  approximate  nonlinear 
strain  displacement  relations  instead  of  the  complicated  exact  ones 
to  describe  T   from  ^T.   This  approach  has  been  considered  by  many 
authors  for  the  nonlinear  analysis  of  thin  shells  [20],  [9],  [2], 
[21]. 


2.2  The  principle  of  virtual  work 

We  consider  the  equilibrium  of  a  thin  shell  structure  in  configuration  T 
subjected  to  body  forces  f  only.   The  internal  stresses  are  described  by 
{  o  }  which  is  a  vector  of  three  components  only,  under  the  assumption  of 
plane  stress  (and  neglecting  the  influence  of  transverse  shear 
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deformation).   The  conditions  of  equilibrium  in   r  leads  to  the  following 
expression  of  the  principle  of  virtual  work: 

y=/<6e>{a}dv-/<6u>{?}dv  =   0  (1) 

v  v 

for  any  {  6  u  }  such  that : 
{  5  u  }  =  {  0  }  on  Su; 
v  is  the  volume  in  T, 
Su  is  the  surface  with  prescribed  displacement, 

{  6  u  }  is  an  arbitrary  virtual  displacement  vector  which  is 
kinematically  admissible, 

{  6  e  }  is  a  virtual  strain  displacement  vector  compatible  with  {  6u  }• 

We  note  that  the  components  of  {  6  u  }  are  defined  with  respect  to  the 
tangent  and  normal  reference  axes  of  the  shell. 
The  3D  virtual  displacement  field  is  of  the  form: 

{  6  u  }  =  {  6  Uj,,  }  +  z  {  6  9  }  (2) 

z  is  the  coordinate  along  the  thickness  h  such  that: 


^  <  z  <  +|  (3) 


{  6  i^  }  are  virtual  displacements  along  two  tangent  directions  x,  y  on 
the  deformed  middle  surface  S  and  along  the  normal  of  S  (Fig.  2).   The 


first  two  components  of  {  5  9  }  are  the  virtual  rotations  around  axes 
tangents  in  S.   The  third  component  of  {  6  0  }  is  zero: 


fS  u^  (6  w, 

{«%}H{v         ;       {66}=   -<6_w,y 


,6 


Co 


Expression  (2)  is  compatible  with  the  so  called  Kirchhoff  hypothesis 
("normal  remains  mormal").  Eqs.  1  to  3  after  integration  through  the 
thickness  give: 


'F=/(<6um>{N}+<6<>{M})d 
s  v 

-  /  (<  6  ^  >  {  fm    }  +  <  6  9  >  {  m  })  d  s  =  0    (4) 
s 

where : 

4 

{  N  }     is  the  three  components  vector  of  membrane  (direct)  forces 

{  M  }     is  the  three  components  vector  of  bending  moments 

{  6  em  }  is  the  three  components  vector  of  virtual  membrane 

strains 

{  6  x  }   is  the  three  components  vector  of  virtual  curvatures 
{  fm  }    are  distributed  forces  along  the  tangent  directions  x,  y 

and  along  the  normal  z 

{?m  )  myy\  (5) 

{  m  }     is  a  two  component  vector  of  distributed  bending  moments  acting 
on  the  shell  surface.   (In  general  these  components  are  zero.) 

We  note  that  expression  4  is  very  general  and  valid  for  any  curved  shell 
surface  where  x  and  y  are  not  necessary  orthogonal  curvilinear 
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Figure  3.   Shell  stress  resultants 
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coordinates.   In  the  case  of  arbitrary  curvilinear  coordinates  all 
expressions  must  be  expressed  in  tensorial  notation  (with  covariant  and 
contravariant  quantities). 

{  6  e^j  }  and  {6k}  are  expressions  in  terras  of  {  6  1%  }  and  of  the 
curvatures  of  the  middle  surface  (  {  6  <  }  is  an  expression  of  the  second 
derivative  of  6  w) . 

The  positive  components  of  {  N  }  and  {  M  }  are  given  for  an  orthogonal 
coordinate  system  on  Figure  3. 

The  stress  resultants  are  related  with  the  stresses  by: 


h  h 

2  2 

{  N  }  =  /  {  a  }  dz  ;  {  M  }  =  /    {  o  }  z  dz  (6) 

-h  -h 

2  2 


The  Euler-Lagrange  expressions  associated  with  the  variational  principle 
(Eq.  4)  are  the  exact  shell  equilibrium  equations  and  the  mechanical 
boundary  conditions  in  terms  of  the  stress  resultants.   These  equations 
are  complicated,  with  coupling  between  {  N  }  and  {  M  }  if  the  shell  is 
described  with  arbitrary  curvilinear  coordinates. 

If  the  shell  is  flat  (or  considered  so)  the  three  equilibrium  equations 
are  the  classical  ones: 


Nx,x  +  Nxy,y  +  fx  =   0 

Nxy,x  +  Ny,y  +  fy  =  0  (7) 

M  +2M  +M  +  T      =    0 

1  lx ,  xx  T   a  Lxy ,  xy   T   *■  y ,  yy  T  L  z        w 


Other  expressions  of  the  principle  of  virtual  work  (Eq.  1  or  7)  can  be 
defined  using  other  reference  configuration  than  T:      the  TLF  involves  °  r 
and  ULF  involves  lT. 
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2.3  The  Total  Lagrangian  Formulation 

2.3.1   General  expressions 

If  we  consider  the  arbitrary  displacement  field  in  Eq.  1  as  the  variation 

■*■ 
of  the  displacement  field  °u  between  °  r  and  T: 

{  6  u  }  =  6  {  0u  }  =  {  6  0u  } 

then  Eq .  1  becomes: 

y=/<60e>    {a}dv-/<    S0u  >    {  T   }  dv  =  0 

v  V     V    {    60u    }  =    {  0    }  on   Su  (8) 

where  the  components  of  {  6  0u    }  can  be  defined  with  respect  to  the 
deformed  (unknown)  coordinates  x,  y,  z  of  V   or  with  respect  to  the 

coordinates  °x,  °y,  °z  in  ° r. 

Eq.  8  can  be  modified  as: 

<F  -   /  <  So  e  >  {  0S  }  d°v  -  /  <  60u  >  {  0T  }  d°v  =  0  (9) 

°v  °v 

where  {  0S  }  are  the  components  of  2n°  Piola  -  Kirchhoff  (P.K.)  stresses 
(tensor  [  0S  ])  and  {  60  e  }  are  the  variation  of  the  Green-Lagrange 
strains  (tensor  [  0e  ]).   [3] 

We  have  the  following  relations: 

{  o?  }  =  oJ  {  T  }  (10) 

[  oS  ]  =  0J  [  0U  p1  [  a  )  [  0u  ]  ~T  (11) 

[  6„e  ]  =  [  SoU  ]T  [  0U  ]  (12) 

J  is  the  Jacobian  of  the  deformation,  i.e.: 


Jo  =  j^  =   det  [  0F  ]  =  det  [  0U  ]  (13) 

where  [  «^  ]  is  the  deformation  gradient  at  a  point  of  the  shell. 

11 


[  0F  ]  can  be  decomposed  as: 

[  oF  ]  =  [  oR  ]  [  oU  ]  (14) 

where  [  0R  ]  corresponds  to  a  pure  rotation  between  a  set  of  coordinates 
in  °T  and  the  deformed  coordinates  in  r  (attached  at  the  same  material 

point).   [  0U  ]  is  the  symmetic  stretch  matrix  for  the  material  point. 
The  Green-Lagrange  (G-L)  strains  are: 

2[  0e  ]  =  [  oF  ]T  [  0F  ]  -  [  I  ]  =  [  0U  ]T  [  U  ]  -  [  I  ]  (15) 

since  [  0R  ]  is  orthogonal. 

The  components  of  [  0  £  ]  are  quadratic  expressions  in  terms  of  the 
component  displacement  of  u  with  respect  to  the  coordinates  x,   y,   z 

o  o     o     o 

of  °T.   They  are  invariant  with  respect  to  rigid  body  motion. 

Eqs.  11,  13  and  15  show  that  under  the  assumption  of  small  strains  we 

have: 

[  oS  ]  =  [  a  ]  (16) 

oJ  «  1  or  d°V  =  dV  (17) 

Eq .  16  means  that  with  the  approximation  of  small  strain  the  2nd  P-K 
stresses  which  are  work  conjugate  to  the  Green-Lagrange  strains, 
corresponds  to  the  "true"  Cauchy  stresses  in  the  deformed  shell.   The  2nd 
P-K  stress  is  therefore  a  material  or  co-rotational  stress. 


This  important  result  is  valid  for  arbitrarily  large  rotations  of  the 
shell  and  will  be  used  for  both  the  TLF  and  ULF  formulation. 

Eqs.  8,  9  and  16  show  also  that: 

{  60e  }  -  {  60e  }  (18) 

12 


Equation  9  is  an  expression  of  terms  of  the  displacements   u  and  of  the 
coordinates  °x,  °y  and  °z  of  °T. 
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2.3.2.   Finite  element  discretization 

We  consider  a  shell  structure  In  its  initial  position  °  r.  This  shell 
will  be  discretized  by  finite  elements  (an  example  using  flat  faceted 
triangular  and  quadrilateral  shell  elements  is  presented  on  Figure  4)  . 

Eq.  9  gives  with  Eqs.  16  and  17: 

*  =  Z   (  /  <  6o  e  >  {  a  }  dve  -  /  <  6cu  >  {  7  }  d  ve  =  0  (19) 

e   ve  ve 

where  dVe  =  d°Ve  represents  the  elementary  volume  on  a  given  element  e. 

If  the  nodal  variables  on  an  element  are  {  0un  }  ,  than  we  can  write: 

{  o£  }  =  [  oB  ]  {  ot^  }  (20) 

{  60e  }  =  [  0B6  ]  {  Soi^  }  (21) 

where  [  0B  ]  and  [  0B6  ]  both  depend  upon  {  0un  } 

Eqs.  19  to  21  give: 

V  =  -  E  <  60un  >  {  0rn  }  =  0  (22) 

e 

with 

{  0rn  }  =  {  0fext  >  "  t  ofint  !  (23) 

{  ofint  }  =  /   [  0B6  ]T  {  a  }  d  ve  (24) 

ve 

/  <  Sou  >  {  T  }  dv  =  <  60un  >  {  0fext  }  (25) 

v 

After  assemblage  of  all  the  elements: 
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Figure  4  A  Shell  discretized  in  Triangular  and  Quadrilateral 
Flat  elements  (from  20) 


15 


*  -  -  <   60Un  >    {  0R   }  =  0 


(26) 


for   all    {    60Un    }  =    {0    }  on   Su 


leading   to: 


{  oR    }  =    {   oFext}  -    {  „Fint    }  =    {  0    } 


int 


(27) 


{  oUn  }  is  the  vector  of  global  nodal  variables  (displacements  between 
°T  and  D. 

{  0R  (0Un)  }  is  the  so  called  residual  global  vector. 

{  o^ext  (  o^n)  J  *s  the   vector  of  the  global  external  forces  that  may  be 
path-deformation  dependant. 

{  oFint  (  oUn)  }  is  the  vector  of  global  internal  forces. 


A  solution  vector  (  0Un  }  is  such  that  {  0R  (  0Un)  }  =  {  0  }  which 
represents  a  system  of  nonlinear  algebraic  equations.   These  equations 

will  be  solved  using  algorithms  and  strategies  based  on  the 

Newton-Raphson  method.   We  need,  therefore,  to  define  a  Jacobian  or 

tangent  stiffness  matrix  [  oK-p  ]  which  results  from  the  assemblage  of 
element  [  ok^  ]  matrices. 

A  symmetic  [  ok^   ]   matrix  is  defined  by  considering: 

«^»f     (<60e>    {    <5   a   }  +  <    6?e  >    {a}-<    S0u>    {   «    })dve        (28) 
ve 

The    first    term   can   in   generai   be   expressed   as: 


II  -  <    60un   >    ([    0k£   ]    +    [    0knJl   ])     {    SoUn    } 


(29) 


The  second  term  is 

I2  -  <  SoU^  >  [  ok0  ]  {  Soi^  } 
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(30) 


The    third   term  is: 

I3  =   <    6oUn  >    ([    0kx  ]     {    5oUn    }  (31) 

So   that  we   have: 

6Y  =   <    Soiin  >    [    okt    ]     {60un    }  (32) 

with 

[   okt   ]   =   [   okA  ]   +   [    o^it  ]   +  [   .kff  ]   -   [   .kx  ]  (33) 

[  0kl   ]    depends  oniy  on  {  01%  }  if  the  material  is  nonlinear. 

[  Oko  ]  is  the  so  called  geometric  stiffness  matrix  of  the  form 

[  .ko  ]  -  /   [  B<|>  ]T  [  N  ]  [  B  <f)  ]  dS  (34) 

se 

where  [  B  <}>  ]  is  constant  in  {  o^  }  and  [  N  ]  is  a  2  by  2  matrix  of 
membrane  forces. 

[  0k^  ]  exists  if  the  loading  is  path  dependant  (the  case  of  hydrostatic 
pressure) . 

In  section  3  the  above  matrix  quantities  are  discussed  for  a  triangular 
flat  faceted  shell  element  within  the  approximation  of  small 
strains  and  moderate  rotations. 
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2.4   The  Updated  Lagrangian  Formulation 

2.4.1   General  expressions 

We  assume  now  that  an  intermediate  configuration  * r  is  obtained  (we  then 
know  all  the  quantities  regarding  geometry  as  well  as  internal  stresses). 

,     1 
(The  internal  stresses  in  l V    {   a  }  are  in  equilibrium  with  the  body 

forces  {  If  }  .  ) 

We  consider  again  Eq .  1  expressing  the  equilibrium  in  Twith: 

{  6  u  }  =  6  {  lU  }  =  {  6ju  }  (35) 

where  {  ^u  }  are  the  displacement  components  from  *■  Y   to  V.      Hence  we 
have : 

i|;  =  /  <  6Le  >  {  a  }  dv  -  /  <  5xu  >  {  f  }  dv  =  0 

V  V 

(36) 
V  {  5]U  }  =  {  0  }  on  Su 

Eq .  8  is  expressed  in  terms  of  2n"  P-K  stresses  with  reference  to  *■ F  and 
variation  of  Green-Lagrange  strain  between  * r  and  T: 

i\>  =   /  <  5ie  >  {  S  }  d*v  -  /  <  5iu  >  {  T  }  dK   =   0  (37) 

1  1  1  1 

with 


{  {¥   }  =  XJ  {  T  }  (38) 

{  !S  }  -  iJ  [  tU  ]"1  [  a  ]  [  iU  J-T  (39) 

{  6lE  }  =  [  5jU  ]T  [  LU  ]  (40) 

!J  =  ^  =  det  [  jF  ]  =  det  [  iU  ]  (41) 
dW 
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[  if   ]  =  [  iR  ]  [  lU  ]  (42) 

2  [  l£  ]  =  [  iF  ]T  [  lF  ]  -  [  I  ]  -  [  XU  ]T  [  XU  ]  -  [  I  ]  (43) 

With  the  approximation  of  small  strains  (section  3.1)  Eq.  36  can  be 
replaced  by: 

$  -  /  <  6je  >    {a}dv-/<    6xu>    {  T   }  dv  =  0 

V  V 

(44) 
V  {  6xu  }  =  {  0  }  on  Su 


with 


{  a  }  =  {  1a   }  +    {  !<J  } 


1 

(45) 
{  f  }  =  {  lf    )  +    (  if  } 

where  {  ^a  }  are  the  incremental  stresses  between  *■  T  and  T  and  {  jf  }  the 
incremental  forces  between  * r  and  T. 

In  Eq .  44  {  S^e  }  are  dependent  upon  the  variations  of  the  displacements 
between  *r  and  r  and  not  upon  the  displacements  {  0u  }. 

In  general  'r  is  a  curved  surface  and  the  exact  expressions  of  [  ^  e  ]  are 

not  simpler  than  [  0e  ].   In  fact  they  are  theorically  identical  when  the 
lower  left  index  o  is  replaced  by  index  1. 
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2.4.2   Finite  Element  discretization 

The  shell  in  configuration  *T  is  discretized  by  finite  elements.   The 
finite  element  matrices  are  of  the  same  general  nature  as  for  the  TLF. 
We  just  have  to  substitute  index  o  with  index  1. 

In  defining  the  tangent  stiffness  matrix  [  ^kt  ]  by  considering  6f  we 
simply  have  to  take  account  of  the  fact  that: 

{  6a  }  =  {  6la    } 
and       {  6f  }  =  {  6{f    } 

In  the  ULF  the  "solution"  means  to  find  the  displacements  and  the 
additional  stresses  between  *r  and  r  that  are  such  that: 

<  1*  >  =  {  lW  >  -  C  lfint  >  =  {  0  }  (46) 
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3.    DESCRIPTION  OF  A  FLAT  TRIANGULAR  SHELL  ELEMENT 

3.1   The  DKT18  shell  element  for  linear  elastic  shells. 

The  nonlinear  analysis  presented  in  this  report  is  based  on  the 
discretization  of  shells  by  flat  triangular  shell  elements  having  three 
nodes  and  the  six  engineering  d.o.f.  per  node  (Fig.  5): 

<  un  >  =  <  Ux  V]_  Wj   RXX   RYj   RZj_ 

U2   V2  W2  RX2   RY2   RZ2  (47) 

U3  V3  W3   RX3  RY3   RZ3  > 

Oj,  vi>  W^   i  =  1,3  are  the  translational  d.o.f  with  respect  to  the 
global  coordinates  axes  X,  Y,  Z. 

RX^,  RY^,  RZ^   i  =  1,3  are  the  rotational  d.o.f.  around  the  global  axis 
X,  Y,  Z. 

The  DKT18  element  results  from  the  superposition  of  the  low  order 
membrane  constant  strain  triangular  element  CST  with  6  d.o.f.  and  of  the 
efficient  bending  triangular  element  DKT  having  9  d.o.f.  (Fig.  6). 

The  linear  stiffness  matrix  of  a  DKT18  shell  element  can  be  expressed 
with  respect  to  the  local  d.o.f.: 

<  "n  >  =  <  ul   vl   wl   9x1   9yl   9zl 

u2  v2  w2   9x2  Qy2      9z2  (48) 

u3  v3  w3   9x3   ey3   ez3  > 
as  [22];  [11]: 

t  ^  ]  =  [  l^n  ]  +  [  kj,  ]  +  [  kez  ]  (49) 
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I  z,  w 

RZ 


RY 


Y,  V 


local  d.o.f.  at  node: 
u,  v,  w,  "'x,  6y,  0z 
global  d.o.f.  at  a  node: 
U,  V,  W,  RX,  RY,  RZ 


x.  u 


Figure  5.   DKT18  Shell  element 
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Figure  6.   CST  and  DKT  plate  elements 
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[  k_  ]  is  the  stiffness  matrix  of  the  CST  element.   It  is  a  simple  matrix 
with  constant  terms.   No  numerical  integration  is  necessary  since  the 

displacement  u  and  v  are  linear. 

[  kjj  ]  is  the  stiffness  matrix  of  the  9  d.o.f.  DKT  plate  bending  element. 
This  element  is  well  documented  in  [6],  [7],  [8]  and  is  obtained  from  the 

the  technique  of  discrete  Kirchhoff  constraints.   This  simple  plate 

bending  element  satisfies  all  convergence  criteria  (like  the  patch-tests) 

and  has  been  found  very  effective  and  reliable  for  thin  plate  bending 

analysis.   It  has  shown  good  behavior  with  respect  to  element 

distortions.   The  stiffness  matrix  of  the  DKT  element  is  obtained  exactly 

(in  linear  analysis)  with  3  numerical  integration  points  in  the  elements 

[19]. 

[  kQz  ]  is  a  fictitious  stiffness  matrix  with  non-zero  components  related 
to  9^1 ,  9z2>  ^Z3  on-ly*   This  matrix  is  necessary  in  order  to  avoid  the 

singularity  of  the  stiffness  matrix  in  the  case  of  copianar  elements. 

The  coefficients  of  this  matrix  should  be  small  enough  so  that  they  do 

not  modify  the  correct  solution  (with  membrane  and  bending  energy  only) 

and  big  enough  to  avoid  numerical  errors.   Two  approaches  are  considered. 

The  first  is  described  in  [  23,  Eq.  13.18  ]  with  a=  10~4  for  our 

computations  on  double  precision  VAX  computers.   The  second  method  is  to 

consider  only  diagonal  coefficients  with  values  a  times  the  minimum  of 

the  diagonal  rotational  coefficients  of  the  bending  stiffness  matrix. 

[10] 

In  the  case  of  symmetrical  material  properties  with  respect  to  the  middle 

surface  of  the  shell  the  stiffness  matrices  [  k^  ]  ,  [  kf  ]  ,  [  k^  ]  are 
not  coupled  so  that  a  large  number  of  coefficients  of  [  k»  ]  are  zeros. 
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The  local  coordinate  x,  y,  z  of  an  element  are  shown  on  Fig.  5.   x  axes 
coincides  with  side  1-2  (origin  in  1)  z  is  normal  to  the  plane  123  (with 
direction  resulting  from  the  cross  products  of  12  x  13).   y  is  such  that 
x,  y,  z  are  orthogonal  and  right-handed.   The  relation  between  the  local 
coordinates  and  the  global  ones  are  [22],  [23]; 


{  x  }  =  [  X  ]  {  X  } 


(50) 


with 


<x>=<x,  y,  z> 
<X>=<X,  Y,  Z> 


(51) 


[  X  ]  is  a  3  by  3  matrix  of  the  direction  cosines  of  x,  y,  z  with  respect 
to  X,  Y,  Z. 

The  element  local  d.o.f .  {  u^  }  are  related  to  the  global  ones  {  un  }  by: 


{  u^  }  =  [  T  ]  {  un  } 


(52) 


with 


[  T  ]  = 


[  X  ] 


0 


[  x  ] 


(53) 


Therefore,  the  stiffness  matrix  of  a  shell  element  in  the  global 
coordinate  system  is: 


[  kjt  ]  -  [  T  ]T  [  kjj  ]  [  T  ] 


If  I  f n  ]  is  a  force  vector  resulting  from  the  discretization  of 
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(54) 


distributed  loads  on  the  elements  (components   fx,  fy  or  fz)  then  the 
corresponding  force  vector  in  the  global,  coordinate  system  is: 

{  f  }  -  [  T  ]T  {  f  }  (55) 

After  the  process  of  assemblage,  modification  due  to  boundary  conditions 
and  solution  of  the  linear  system  we  can  obtain  the  strains  and  the 
stresses  at  any  point  in  the  element: 


{  H   >  =  I  Bm  1  (  u™  }  +  z  [  Bb  ]  {  u£  }  (56) 


~n  u  ~n 


where 


<  uJJ  >  =  <  Ul  vj   u2  v2  u3  v3  > 

<  H*  >  =  <  wl   9xl   8yl  w2   9x2   9y2  w3   ^3   ^3  >  <57> 


[  T^  ]  is  constant  and  [  B^  ]  (the  linear  strain  operator  of  the  DKT 
plate  element)  is  linear  in  x,  y. 

In  the  absence  of  coupling  between  membrane  and  bending  effects  we  have: 

{  N  }  =  [  Dm  ]  [  Bm  ]  {  u£  } 

(58) 
{  M  }  =  [  Db  ]  [  Bb  ]  {  ub  } 

and  therefore : 

{  a  }  =  [  1^  ]  [  Bm  ]  {  u£  }  +  z    [    %    }     [    Bb  ]  {  u£  }  (59) 

[  Dm  ]  and  [  D^  ]  are  3  by  3  membrane  and  bending  material  matrices.   In 
the  ususal  case  of  plane  stress  isotropic  material: 
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Eh 


[  Dm  1  =  l-v2 


1   v  0 

v  1   0 

0   0  j^ 
2  -J 


12 


(60) 


where  E  and  v  are  Young's  modulus  and  Poisson's  ratio. 

Eq.  58  shows  that  in  general  the  membrane  forces  are  constant  and  the 
bending  moments  vary  linearly. 

We  have  computed  the  stress  resultants,  the  principal  stresses  and  the 
Von  Mises  equivalent  stress  on  the  outer  faces  of  the  shell  at  the 
maximum  of  7  points  per  element  (centroid,  integration  points,  corner 
nodes).   The  corner  node  values  are  discontinuous  but  can  provide  useful 
information  with  respect  to  node  location  and  with  respect  to  precision 
in  the  results. 

The  simple  DKT18  shell  elements  has  been  used  extensively  for  linear 
analysis  of  shells  and  is  implemented  in  several  computer  codes  working 
on  mini  and  micro-computers.   The  main  disadvantage  of  the  element  are 
the  CST  element  as  membrane  element  and  sometimes  the  non-energy 
associated  0£  d.o.f.'s. 
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3.2  The  element  matrices  for  TLF  nonlinear  analysis 

Our  TLF  combined  with  the  use  of  the  flat  triangular  DKT18  shell  element 
is  based  on  the  following  definition  of  the  G-L  strain: 

{  oe  }  -  {  h  >  +  {  %*  >  (61) 

with 

u,x  +  z  9y,x 

{  ei   ^   =  {  v'y  "  z  ^'y        (  ^62^ 

u»y  +  v»x    "  z  ^'x  "  Vy 

1  2 

2  w>x 

{    ^i    }  =  (I    w,|  }  (63) 

w,x  w,y 

where  the  lower  left  index  o  has  been  omitted  everywhere  for 

simplification,  i.e.,u,  v,  w,  9jj  and  9y  are  displacements  and  rotations 

with  respect  to  axes  x,  y,  z  of  the  undeformed  shell  element. 

{  ££  }  is  the  vector  of  linear  strains  which  leads  to  the  linear 
stiffness  matrices  [  k^  ]  and  [  k^  ]  presented  in  section  3.1. 

{  e^  }  involves  only  derivative  of  w  with  respect  to  x  and  y.   This 
nonlinear  part  of  G-L  strain  is  associated  with  the  so-called  Von  Karman 
plate  theory,  that  is  this  expression,  will  be  valid  only  for  large 
displacements  and  moderate  rotations.   Therefore,  the  TLF  discussed  here 
is  valid  with  the  above  assumption  and  a  flat  triangular  discretization 
of  the  shell  in  its  initial  configuration. 

Eq .  61  gives  : 

{  60e  }  =  [  B6  ]  {  6  uj,  }  (64) 
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where  <  61^  >  =  <5  <  "n  > 
with  <  Un  >  given  in  Eq.  48 
[  B6  ]  -  [  B*  ]  +  [  B^  ] 


(65) 


[  B^  ]  is  the  linear  operator  which  leads  to  [  k™  ]  and  [  k^  ] 


t  BfU  1  " 


w 


w,x   <  N  x  > 


w 


w,y  <  N  y  > 

» 

w,y  <  NWX  >  +  wlx  <  NW  > 


with 


(66) 


w  =  <  N"  >  {  u^  }  (67) 

<  N*7  >  has  nine  non  zero  components  which  are  associated  with  a  cubic 
Hermite  interpolation  function  for  w.   This  9  term  interpolation  function 
is  chosen  so  that  it  is  invariant  with  respect  to  local  coordinate  x,  y. 
This  incomplete  cubic  interpolation  is  given  in  [23,  Eq .  10.29]. 
The  internal  forces  Eq .  24  for  the  element  in  the  global  coordinate 
system  is  then  defined  as  : 


{  ofint  }  =  [  oT  ]T  /   [  B6  ]T  {  a  }  dve 


(68) 


with  [  0T  ]  -  [  T  ]  given  in  Eq .  53  and  [  B6 
the  case  of  elastic  behavior: 

{    a    }  =    [    D   ]    ([    BZ   ]    +1   [    Bni  ])     {u^    } 


]    given   in   65  and   66.      For 


(69) 


{  0fext  }  ^ *  25  depends  on  the  loading.   If  the  loads  are  not  path 
dependant  {  fext  }  is  constant,  if  not  {  fext  }  is  a  function  of 

{  Uj,  }.   The  interpolation  functions  for  the  evaluation  of  the  equivalent 
forces  (Eq.  25)  are  considered  linear  for  u,  v,  w. 
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The  following  expression  can  be  considered  in  the  case  of  uniform  normal 
pressure  of  intensity  p  with  respect  to  the  deformed  middle  surface: 


6Wp  =  p  /  (  -  w,x  6u  -  w,y  6v  +  6w  )  dxdy 
se 

=  <  «an  >  <  f4xt  }  "  <  6un  >  <  fJxt  } 


[  -falA  1  -  /   [  B£  ]T  [  D  ]  [  Bn£  ]  +  [  BnJl  ]T  [  D  ]  [  Bji 
ve 

+  [  BqA  ]T  [  D  ]  [  Bnil  ]  dVe 


Moreover: 
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(70) 


with 


{  ofPxt  }  -  [  T  ]T  {  fP^  }  (71) 


The  tangent  stiffness  [  0kt  ]  as  defined  in  Eq.  33  is  such  that: 

[  okt  ]  =  [  T  ]T  [  okj.  ]  [  T  ]  (72) 

where 

I  okt  1  =  [  H.Z  1  +  I  °ki£  1  +  [  oka  ]  (73) 

[  k^  ]  is  given  in  Eq.  49. 
For  elastic  material: 

{  6a   }  =  [  D  ]  {  6e  }  (74) 

Therefore: 


(75) 


f  °£na  ]  =  /   [  B.  ]T  [  N  ]  [  B.  ]  dxdy  (76) 


where 


[  N  ]  - 


N. 


xy 


Nxy   Ny 


(77) 


[  B«|)  ]  = 


w 
U<  N,y  >J 


(78) 


with  <  Nw  >  the  9  term  cubic  function  Eq.  67. 

The  [  k^  ]  matrix  is  not  necessary  for  our  moderate  rotation  TLF. 

Matrices  [  okj^  ]  and  [  ok0  ]  are  evaiuated  with  no  negiecting  terms  in 
the  tangent  stiffness  matrix  using  3  numerical  integration  points. 

The  TLF  as  presented  above  will  give  the  nonlinear  solution  of  arbitrary 
shell  structures  within  the  approximations  considered.   That  is,  the 
converging  solution  with  mesh  refinement  will  always  be  restricted  to  the 
moderate  rotation  assumptions. 
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3.3  The  element  matrices  for  ULF  nonlinear  analysis 

In  our  Updated  Lagrangian  Formulation  using  the  DKT18  shell  elements  the 
intermediate  configurations  ^T   are  not  the  exact  configurations  of  the 
shell.   These  approximate  configurations  will  result  from: 

-  the  discretization  of  the  shell  using  flat  triangular  elements 

-  the  assumption  of  moderate  rotations  between  two  configurations  (like 
between  *T  and  T) . 

Hence  the  configuration  T'is  obtained  from  * T  by  making  the  same 
assumptions  and  the  same  type  of  computations  as  between  °  r  and  r  in  the 
TLF  procedure. 

So  we  assume  that  the  current  coordinates  in  *■  V   are  known.   They  result 
from: 

{  U    }  =    i   0x  }  +  {  J  u  }  (79) 

(The  curvatures  in  1 r  are  neglected  as  they  are  in  °  r) .   We  also  assume 
that  the  Cauchy  "true"  stresses  {   a  }  are  obtained  (and  stored  at  the 
integration  points  of  the  triangular  elements).   These  stresses  are  in 

equilibrium  with  the  surface  forces  0fx  ,  0fy  and  0fz« 

The  necessary  information  to  obtain  r  are  the  residual  vectors  and  the 
tangent  matrices  of  each  element.   These  quantities  are  obtained  in  a 
similar  manner  as  in  TLF. 

The  internal  forces  {  ifint  )  are  defined  as: 

<  lfint  >  =  [  tT  ]T  /   [  jBd  ]T  {  a  }  dV  (80) 

ire 
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[  iT  ]  is  the  matrix  of  direction  cosines  between  the  iocai  axes  *x,  *y, 
*z  in  *T  and  the  global  axes  X,  Y,  Z. 

[  jB6  ]  is  similar  to  Eqs.  65  and  66  where  w  is  the  displacement  in  the 
*z  direction  between  *r  and  r. 

{  a  }  is  defined  as: 

{  a  }  =  [  jo  }  +  {  xa  }  (81) 

and 

{  i<3    }  =  [  D  ]  ([  Bz   ]    +1  [  1Hnl  ]     {  xun  }  (82) 

{  ^un  }  are  the  nodal  d.o.f.  in  the  global  coordinate  system  and  refer  to 
displacements  and  rotations  between  *■  T   and  T. 

The  external  forces  {  ifext  }  include  the  forces  from  "r  to  I*. 

The  tangent  stiffness  matrix  is  kept  complete  and  is  given  by: 

[  l^t  1  »  f  lT  ]T  [  ikj-  ]  [  iT  ]  (83) 

where  [  ^kt  ]  is  similar  to  [  0kj-  ]  . 

We  again  note  that  the  geometric  stiffness  matrix  [  \k_0  ]  contains  the 
influence  of  the  stresses  from  °  T  to  r,  and  that  [  i^^  ]  is  nonlinear 
in  terms  of  ^w  (from  1 r  to  D . 

We  note  that  the  follower  forces  are  easily  taken  into  account  in  the  ULF 
since  the  coordinates  (and  therefore  the  element  orientations)  are 
updated  after  each  new  known  configuration. 

The  performance  of  the  ULF  in  computing  with  precision  the  nonlinear 
response  of  shells  with  large  rotations  will  not  only  depend  on  the 
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number  of  elements  but  also  on  the  number  of  steps  (or  configurations  *  r) 
between  °Y   and  the  unknown  configuration  r,  because  of  the  assumption  of 
moderate  rotation  between  two  configurations. 
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4.    ON  THE  AUTOMATIC  SOLUTION  FOR  PRE  AND  POST  BUCKLING 

4.1   Solving  the  nonlinear  equations. 

The  solution  in  both  TLF  and  ULF  must  satisfy  a  set  of  simultaneous 
nonlinear  algebraic  equations  as  given  by  Eq.  27  or  46  of  the  form: 

{  R  (U,X)  }  =  X  {  Fext  (U)  }  -  {  Fint  }  =  {  0  }  (84) 

where  the  number  of  equations  n  is  equal  to  the  total  active  d.o.f.  of  the 
discretized  problem,  ({  U  }  stands  for  these  active  d.o.f).   X  is  a  load 
parameter  (we  consider  only  one  variable  loading).   {  F^nt  }  is  always  a 
function  of  {  U  }  and  {  FexC  }  is  so  only  if  the  loading  is  path 
dependent . 

In  the  nonlinear  analysis  of  shell  structures,  the  "load-displacement" 
curves  can  exhibit  all  kinds  of  forms  depending  on  the  problem  (geometry, 
loading,  boundary  conditions,  material  properties).   [1],  [2]... 

In  this  report  we  consider  only  elastic  behavior  of  shells  with  large 
displacements  and  large  rotations,  pre  and  post-buckling  with  multiple 
limit  points,  snap-through  and  snap-back  behavior.   The  problem  of 
determining  Euler  bifurcation  loads  by  solving  linear  eigenvalue  problems 
is  not  considered  although  the  basic  ingredients  (stiffness  matrices, 
geometric  stiffness  matrices,  and  eigenvalue  equation  solvers  are 
available).   Bifurcation  loads  can,  however,  be  obtained  after  the 
introduction  of  a  perturbing  parameter  that  discloses  the  bifurcation 
mode. 

The  complete  determination  of  the  load  displacement  curves  can  be 
performed  using  different  strategies  all  based  on  a  number  of  iterative 
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methods.   The  problem  is  to  obtain: 

-  the  n  components  of  {  U  }  for  a  given  X  or 

-  the  n  components  of  {  U  }  and  X  with  one  constraint  equation. 

The  over-all  behavior  of  the  load-displacement  curves  can  be 
characterized  by  the  so-called  current  stiffness  parameter  [24]  Sp. 

One  simple  definition  in  the  case  of  constant  loading  is: 

-  ,  AXP  <  AU0  >  {  F  } 

SP  "  AX^  <  AU^  >  {  F  }  (85) 

p  p 

{  AUjj  }  is  the  linear  solution  for  AX^.   (a  U   }  and  AX  are  the 

increments  of  displacements  and  of  load  at  step  p.   Sp  is  a  useful 

parameter  in  an  automatic  determination  of  the  complete  load  displacement 

curves . 

Three  strategies  have  been  implemented  and  used  to  solve  various 

nonlinear  shell  problems.   The  first  is  the  load  control  strategy 

(prescribed  X),  the  second  is  the  one-displacement  control  strategy  (one 

prescribed  component  of  {  U  }) ,  the  third  is  the  arc-length  strategy 

involving  all  d.o.f.. 

The  three  strategies  are  using  the  Newton-Raphson  method  to  obtain  the 
incremental  solutions  and  are  as  automatic  as  possible  within  their  own 
limitations . 
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4.2   The  Newton-Raphson  method  with  prescribed  forces 


The  algorithm  is  the  following  [  ],  (Fig.  7a) 


step  p   :   \P,  {  UP  }   known  solution 
step  p+1:   XP+1  =  ~X     {  U1  }  =  {  UP  } 

iterations:   i  =  1  to  NITER 

{  R1  }  =  I  {  Fext  (U1)  }  -  {  Flnt  (U1)  } 

[  K^  ]  {  AU  }  =  {  Ri  } 

{  Ui+1  }  =  {  U1  }  +  {  AU  } 
TEST  convergence 
step  p  +  2  ... 


(86) 


where  [  K_  ]  is  the  global  tangent  stiffness  matrix  which  is  computed  at 
each  iteration  in  the  full  Newton-Raphson  (N-R)  method  or  computed  only 

at  the  beginning  of  the  iteration  process  in  the  modified  N-R  method. 
The  test  of  convergence  is: 


TEST  <  EPSILON 


(87) 


with  TEST  defined  as: 
AU  II 


TEST1  = 


or 


U1 


(88) 


TEST2  = 


AU 


II  Ui"  UP  || 


(89) 


where  ||  U    =  (<  U  >  {  U  })  ^2  is  the  Euclidean  norm  of  the  total 

displacement  vector.   We  have  usually  considered  EPSILON  =  10--5  if  TEST1 

is  used  and  10~2  if  TEST2  is  used.   The  TEST2  is  motivated  so  that  all 
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Figure  7   Nonlinear  solution  strategies 
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steps  (particularly  in  the  ULF)  will  have  the  same  convergence 
requirements. 

\   is  fixed  by  the  user  or  determined  automatically.   If  no  convergence 
has  occurred  (i.e.,  i>  NITER  and  TEST  >  EPSILON)  than  the  given  ~X  or  ~AA  is 
cut  by  two  automatically  until  convergence  is  reached.   A\j  =  XP+^  -  )P 
can  also  be  modified  depending  on  the  convergence  rate  and  according  to: 

"^p  -  aVi  YFT  a  (90) 

where  AAp_i  is  the  increment  of  loading  at  the  previous  solution  step. 
Id  is  a  number  of  required  iterations  and  In-\    is  the  number  of 
iterations  at  step  p-1.   a  is  a  number  defined  by  the  user.   We  have 
considered  a  =  1  and  0.5  with  Id  =  4  in  the  full  N-R  and  6  in  the 
modified  N-R. 

The  above  strategy  has  been  found  effective  to  obtain  the  load  deflection 
curves  automatically  up  to  the  first  limit  point,  giving  the  buckling 
load,  (and  therefore  the  complete  nonlinear  response  if  there  is  no 
limit  point). 
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4.3  The  N-R  method  with  a  prescribed  displacement  component. 

This  algorithm  has  been  found  very  efficient  to  obtain  the  post-buckling 
response  when  a  particular  component  of  the  displacement  vector  stiil 
increases  after  the  limit  load.   [12],  [13],  [25] 

The  algorithm  is  the  following  (Fig.  7b): 

step  p    :   XP,  {  UP  }   known  solution 
step  p+1  :   X1  =  XP  ;     {  U1  }  =  {  UP  } 
r- iterations:   i  =  1  to  NITER 


{  R1  }  =  X1  {  F  \  -  {  F.\   } 
ext      int 

[  K^  ]  (  {  AUR  }   {  AUF  })  =  (  {  Ri  }   {  Fe^  }) 


(91) 


{  U 


i+1 


}  =  {  U1  }  +  {  AUR  }  +  AX  {  AUF  } 


X 


i+1  .  ii 


=  X1  +  AX 


where  AX  is  such  that 

(AUR)q  +  AX  (AUF)q  =  AUq 
'  TEST  convergence 


step  p  +  2  ... 

where  AUq  is  a  prescribed  displacement  increment,  (  AUR)q  and  (  AUp)q  are 
the  qth  component  of  vectors  {  AUR  }  and  {  AUF  }.   It  is  also  possible  to 
use  the  modified  N-R  method. 

For  the  same  problem  and  convergence  test  the  above  algorithm  leads  in 

general  to  a  faster  rate  of  convergence  compared  to  the  prescribed 
loading.   This  is  due  to  the  modification  of  [  1C,  ]  after  the  first 
iteration.    However,  two  load  vectors  are  considered  at  each  iteration. 


As  in  the  previous  algorithm,  AUq  can  be  automatically  adjusted  if 
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convergence  doesn't  occur  within  the  limitation  given  by  the  user.   This 
algorithm  is  very  efficient  in  many  snap-through  situations  and  works 
until  a  limit  point  in  displacement  (snap-back)  occurs. 

The  above  algorithm  is  a  particular  case  of  the  arc-length  algorithm  as 
discussed  below. 
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4.4   The  N-R  and  arc-length  control  method. 

The  so-called  arc-length  or  modified  arc-length  method  has  received  a 
great  attention  in  the  last  five  years  [14],  [15],  [16],  [17],  [18], 
[26]. 

The  algorithm  is  similar  to  the  previous  one  (displacement  control).   It 
is  only  different  in  the  evaluation  of  AX  at  the  first  and  subsequent 
iterations.   (Fig.  7c): 


step  p    :   AP,  {  UP  }   known  solution 
step  p+1  :   X1  =  XP  ;     {  U1  }  =  {  UP  } 

r  iterations:   i  =  1  to  NITER 

{  R1  }  =  X1  {  F  *   -  {  F  *   } 
ext      mt 

[  K^  ]  (  {  AUR  }   {  AUF  })  =  {  Ri  }   {  F 


ext 


}) 


ri+1 


{  U1+1  }  =  {  U1  }  +  {  AUR  }  +  AX  {  tfJF  } 


Xi+1  =  X1  +  AX 


where  AX  is  such  that 


<  Ui+1  -  UP  >  {  Ui+1  -  UP  }  =  (ASp)2 
and 

<  U1  -  UP  >  {  U1+1  -  UP  }  >  0   i>l 

<  UP  -  UP-1  >  {  U2  -  UP  }  >  0   1=1 
L  TEST  convergence 

step  p  +  2  ... 


(92) 


(93) 


(94  a,b) 


Eq .  93  is  a  quadratic  equation  in  AX  that  can  be  written  as 


a   AXZ  +  b   AX  +  c   =   0 


(95) 


with 
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a  =  <  AUF  >  {  AUF  } 

b  =  2  <  AUF  >  {  V  } 

c  =  <  V  >  {  V  }  -  (  ASp)2 

<V>  =  <  AUR  >  +  <  U1  -  UP  > 

If  no  real  root  of  Eq.  94  exists, the  arc-length  ASp  must  be  reduced.  The 
choice  of  the  real  root  is  such  that  Eq.  94a  or  b  is  satisfied. 

One  should  mention  here  that  other  definitions  of  the  arc-length  (Eq.  93) 
can  be  made,  but  the  above  relation  has  been  found  effective  to  solve  our 
examples . 

As  in  the  two  other  strategies,  it  is  possible  to  adjust  automatically 
the  value  of  ASp  between  two  steps: 

-  no  convergence   ASp  =  0.5   ASp_i 

Id 


-  if  convergence   ASp  =  ASp_^ 


Vi 

with  Id,  Ip-i>  a  as  discussed  in  section  4.2. 


(96) 


If  the  arc-length  strategy  is  used,  we  first  start  the  problem  using  one 
of  the  two  previous  methods  (load  or  displacement  method).   Then  we 
compute  AS2  using  Eq .  93  to  obtain  the  solution  at  step  2  using  the 
arc-length  algorithm. 
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5.    NUMERICAL  RESULTS 

5.1   Comments  on  the  computational  procedure 

A  FORTRAN  program  has  been  written  partly  for  this  research.   The  basic 
routines  of  the  finite  element  method  are  those  documented  in  [19].   We 
have  made  extensive  modifications  in  the  nonlinear  block  to  implement  our 
methods  and  strategies.   We  have  also  written  the  routines  dealing  with 
the  triangular  shell  element. 

The  examples  discussed  below  are  solved  using  a  VAX  11/780  or  an 
APOLLO/ DN300. 

Only  simple  problems  have  been  solved  and  presented  in  this  report.   They 
involve  a  limited  number  of  d.o.f.  (about  200).   These  examples  are 
chosen  in  order  to  show  the  various  possibilities  of  the  present 
formulation  to  deal  with  pre  and  post-buckling  and  large  rotations  of 
arbitrary  shells. 
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5.2  Nonlinear  response  of  a  3D  truss  structure. 

This  example  taken  from  [27]  is  chosen  to  show  the  possibility  of  the 
implemented  numerical  methods  to  deal  with  the  automatic  computation  of 
very  nonlinear  problems.   The  3D  structure  shown  in  Fig.  8  is  made  of  24 
truss  elements  having  2  nodes  and  three  d.o.f.  per  node  (the  3  displacement 
components  u,  v,  w) .   The  structure  is  fixed  at  the  base  and  subjected  to 
a  point  load  at  the  center  1  (in  the  u  direction  on  Fig.  8).   There  are 
21  active  d.o.f. 

A  large  number  of  runs  have  been  performed  with  various  parameters  such 
as : 

-  TLF  or  ULF  options 

-  arc-length  with  or  without  adjusting  ASp 

-  full  N-R  or  modified  N-R 

-  load  value  at  the  first  step 

-  influence  of  the  TEST  of  convergence  in  the  numerical  process 

Some  load-displacement  curves  are  given  on  Figs.  9  and  10,  where  u  is 
the  displacement  under  the  load  P  and  v  is  the  displacement  in  direction 
y  at  node  2.   Figure  9  is  obtained  with  the  full  N-R  and  the  automatic 
arc-length  method  with  no  modification  of  arc-length  (Eq.  88  is  used  with 
e  =  10~3).   The  first  nonlinear  solution  is  obtained  with  a  prescribed 
value  of  P  =  1CF  P/EA  equal  to  2.   Then  the  arc-length  ASD  is  computed 


using  Eq.  93,  i.e.  ASD  =  <  U1  >  {  U1  }.   All  symbols  in  Figure  9  coincide 


'P 

=     /     III      S       t     TT 1       \  All       clmKn  I  c      -in      T 

with  equilibrium  solutions  obtained  automatically.   There  are  8  limit 
points  for  the  range  of  load  and  displacements  considered  and  these 
curves  correspond  to  the  primary  solution  (with  full  symmetry).   Of 
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Figure  8   Three  dimensional  truss  structure 
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course  there  are  other  solutions  which  are  associated  with  bifurcations. 
Figure  10  presents  the  solution  in  the  case  of  automatic  modification  of 
the  arc-length  using  Eq.  96  where  a  =  0.5  and  Id  =  4.  All  points  on 
Figure  10  correspond  to  solutions.   This  figure  shows  the  robustness  of 
the  automatic  computation  algorithm  when  full  N-R  and  arc-length  methods 
are  combined.   Our  results  coincide  with  those  presented  in  a  recent 
paper  [27] . 
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5.3   Snap  through  the  snap  back  of  the  cylindrical  shell  (CTEX4) 

The  problem  presented  in  Figure  11  has  been  widely  used  in  the  literature 
to  compare  the  performance  of  various  nonlinear  formulations,  finite 
element  models  and  nonlinear  solution  strategies.   [18],  [20],  [28], 
[16]. 

A  simple  mesh  of  DKT18  elements  was  considered  (48  elements,  210  dof 
before  elimination  of  the  imposed  variables).   The  straight  edges  are 
hinged  and  the  curved  edges  are  free. 

In  this  problem  we  have  studied  the  influence  of  the  formulation  (TLF 
verse  ULF) ,  the  influence  of  the  arc-length  strategy  on  the  solution. 
The  solutions  have  been  obtained  using  the  full  N-R  or  the  modified  N-R 
method . 

The  influence  of  the  formulation  can  be  seen  on  Figures  12  and  13  where 
curves  relate  the  load  versus  the  normal  under  the  load  or  at  the  free 
edge.   (The  results  have  been  obtained  with  full  N-R,  with  constant 
arc-length  steps).   The  first  solution  was  obtained  for  a  prescribed  load 
P  =  lkN.   Figures  12b  and  13b  are  in  good  agreement  with  the  "reference 
solution"  as  given  by  several  authors.   The  TLF  gives  a  higher  value  for 
the  buckling  load  and  doesn't  reproduce  the  snap-back  behavior  for  the 
central  displacement  for  the  small  mesh  considered.   However  it  is 
expected  that  the  correct  answer  will  be  obtained  with  the  refinement  of 
the  mesh  since  this  problem  doesn't  involve  very  large  rotations. 


47 


The  curves  on  Figure  14  are  obtained  using  the  modified  N-R  method,  the 
ULF  and  a  variable  arc-length  increment  (with  a=  0.5  and  Id  =  4).  The 
results  are  the  correct  ones  and  are  obtained  very  efficiently  in  terms 
of  computer  time. 
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Figure   11.      Cylindrical    shell  with   free   curved   edges   and   hinged 
straight   edges    (CTEX4) 
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Figure  14   P  versus  W  and  W.  for  CTEX4 
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5.4   Spherical  cap 

The  problem  presented  on  Figure  15  has  also  been  considered  by  several 
researchers  since  1969  [29],  [20],  [30],  [31],  [32].   The  four  edges  are 
hinged  and  a  load  is  applied  at  the  center.   A  uniform  mesh  of  5  by  5 
elements  (50  elements,  216  d.o.f.  before  the  elimination  of  the 
prescribed  d.o.f.)  has  been  considered. 

The  central  deflection  versus  the  load  is  given  on  Figure  16  for  both  TLF 
and  ULF  and  using  the  imposed  displacement  and  the  full  N-R  methods  (Eq. 
91).   The  value  AUq  is  constant  and  equal  to  0.2  h.   The  ULF  gives  good 
results,  in  agreement  with  the  results  given  by  other  authors.   The  TLF 
leads  to  a  slightly  higher  buckling  load  and  doesn't  represent  the 
unstable  branch  properly.   A  finer  mesh  would  result  in  better  results 
using  the  TLF  here  since  the  displacements  are  only  two  times  the 
thickness.   The  same  type  of  results  have  been  obtained  using  the 
arc-length  algorithm  with  or  without  variable  ASp. 
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Figure  16.   Spherical  cap.   Load  versus  central  deflection 
for  TLF  and  ULF 
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5.5   Cylindrical  shell  with  clamped  curved  edges  (CTEX1) 

The  membrane  stiffening  behavior  of  a  cylindrical  panel  subjected  to  a 
central  load  with  straight  free  edges  and  clamped  curved  edges,  as 
described  in  Figure  17,  has  been  studied  with  a  4  by  6  mesh  for  a  quarter 
of  the  shell. 

The  central  displacement  versus  the  load  is  presented  in  Figure  18.   The 
numerical  results  presented  are  obtained  using  load  incrementation 
(constant  increments  AX  =  1  lb)  and  full  N-R,  for  both  TLF  and  IJLF.   With 
the  TLF,  the  results  are  quite  far  from  those  of  [33]  using  sophisticated 
cubic  Lagrangian  isoparametric  elements  (a  mesh  of  4  x  6  elements  leading 
to  1200  d.o.f.  was  considered  in  [33]).   The  results  using  TLF  elements 
cannot  be  improved  by  reducing  the  load  increments  since  convergence  has 
occurred  and  there  is  no  influence  of  the  load  steps  on  the  converging 
results.   These  results  can,  however,  be  improved  by  using  finer  meshes. 
The  results  with  TLF  are  not  good  because  we  have  fairly  large 
displacements  (up  to  10  times  the  thickness). 

For  the  same  number  of  load  steps,  the  ULF  formulation  gives  better 
results  than  the  TLF.   This  is  due  to  the  effect  of  large  displacements 
and  moderately  large  rotations.   Improved  results  can  be  obtained  with 
the  ULF  if  the  load  steps  are  reduced.   The  correct  results  would  be 
obtained  if  both  the  number  of  elements  and  the  number  of  steps  are 
increased. 
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Figure  17.   Cylindrical  shell  with  straight  edges  free 
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5.6  Far  post-buckling  of  a  cylindrical  shell  with  hinged  curved  edges 
(CTEX2) 

The  shell  structure  presented  in  Figure  17  is  again  considered  but  with 
hinged  conditions  on  the  curved  edges  (these  edges  are  not  restrained  in 
the  axial  direction). 

In  this  case  the  behavior  is  different  from  the  clamped  case  since  the 
displacements  are  much  larger  for  the  same  load  and  snap  through  occurs 
for  a  load  of  2.3  lbs. 

The  same  type  of  analysis  as  in  the  clamped  case  has  been  performed  for 
0  <  W/h  <  14  using  now  the  prescribed  displacement  algorithm  for  various 
AUq  =  Awq  and  again  the  full  N-R  method.   The  results  are  presented  in 
Figure  19  for  the  TLF  and  the  ULF  (7  and  14  steps).   The  inability  of  the 
TLF  to  find  a  limit  load  is  clearly  shown  in  the  figure.   Again  better 
results  would  need  much  more  elements.   The  ULF  leads  to  good  results 
with  a  limit  load  20%  higher  than  the  reference  value  taken  from  [33] 
with  28  steps.   A  more  accurate  value  would  need  more  d.o.f.'s. 

The  response  of  the  shell  for  very  large  displacements  and  rotations  has 

also  been  studied  for  the  same  mesh.   Figure  20  shows  P  versus  the 

central  displacement  up  to  150  times  the  thickness  (and  1/4  of  the 

length)  using  70  steps,  the  ULF,  full  N-R  and  the  displacement  control 

algorithm.  The  slight  oscillations  observed  on  Figure  20  are  due  to  the 

fact  that  in  the  regions  considered  the  overall  tangent  matrices  are  very 

ill-conditioned.   The  influence  of  the  type  of  representation  of  [  k^  ] 
in  these  regions  is  important.   If  the  displacement  increments  are 

reduced,  these  oscillations  disappear  and  the  behavior  is  slightly 
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different  as  shown  in  Figure  21  where  140  steps  are  considered.   One  can 
observe  on  this  figure  the  low  post-buckling  minimum  and  a  second 
snap-through  for  W/j  =  90  which  corresponds  to  a  local  deformation  of  the 
cross-section  near  the  boundary.   The  results  presented  in  Figure  21  are 
in  good  agreement  with  those  reported  in  [33] ,  the  most  important 
difference  being  in  the  evaluation  of  the  first  buckling  load  as  shown  on 
Figure  20. 

The  above  results  have  also  been  obtained  using  the  automatic  variable 
arc-length  algorithm  combined  with  the  full  N-R  method.   This  example  has 
clearly  shown  the  capability  of  the  formulation  to  deal  with  large 
rotations . 
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6.    CONCLUSIONS 

The  numerical  results  presented  in  this  report,  have  shown  the  ability  of 
the  DKT18  elements  combined  with  the  ULF  to  give  accurate  and  efficient 
answers  to  very  different  types  of  nonlinear  shell  problems  including 
snap-through,  snap-back  situations  and  lar^e  rotations  post  buckling. 
The  three  automatic  strategies  have  been  tested  and  it  is  found  that  the 
full  N-R  method  combined  with  the  arc-length  method  are  very  reliable  and 
powerful  to  deal  with  all  kinds  of  nonlinear  situations.   The  overall 
package  can  solve  moderately  large  problems  on  mini  and  micro-computers. 
This  package  has  several  capabilities  since  it  has  retained  different 
aspects  such  as  TLF,  ULF,  full  or  modified  N-R,  automatic  constant  or 
variable  load  or  displacement  or  arc-length  increments,...   The  modules 
dealing  with  the  nonlinear  procedures  and  the  DKT18  shell  elements  can 
be  adapted  to  other  finite  element  codes  having  a  similar  structure  than 
MEF  [19].   Also, the  procedures  are  independent  of  the  DKT18  shell 
element i  therefore,  the  library  of  elements  can  be  enriched  in  the 
future.   Elasto-plastic  behavior  can  also  be  included  if  small  strains 
are  assumed. 
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